# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
#       Patriotims!
#
#       FIRST VERSION  June       3, 2022
#       THIS VERSION   June       9, 2022
#       LAST RUN       June       9, 2022
#
#       LAST REVISOR   BC
#
#       Log of revisions:
#
#       This creates all the maps we need for the project.
#       In particular:
#          - a CountyND with the boundaries for the ND counties
#          - the intersections between the CountyND & all other maps
#       Starting counties are 1930. Relative to the original, the new shapefile:
#               - Has 1 observation for each ND county
#               - AREA (in m2)
#               - CENTROID COORDINATES (with WGS 1984)
#
#                          PLAN OF THE PROCEDURE
#           1.      Build the New Deal shapefile
#           2.      Compute intersections
#                 a. With New Deal counties
#                 b. With 1940 counties
#                 c. With Climatic divisions
#                 d. With Electoral districts
#           3.      Erase junk
#
#        The file was written and run on Python 2.7.9 and ArcMap 10.6.1
#
# ---------------------------------------------------------------------------

# Import arcpy module
print "I am importing the libraries."
import arcpy
import arcgisscripting
import os
print(arcpy.ListVersions())
#gp = arcgisscripting.create()
#gp.OverWriteOutput = 1

# ---------------------------------------------------------------------------
#           1.      Build the New Deal shapefile
# ---------------------------------------------------------------------------
# Local variables:
print "              Part 1: create the shapefile of the New Deal counties"
print "I am defining the variables."

cwd  = os.getcwd()
path = cwd.replace("code-build","")

try:
    os. mkdir(path + "/tmp")
except:
    pass

dbf     = r"" + path + "/rawdata/maps/Bridge.dbf"
shp0    = r"" + path + "/rawdata/NHGIS/US_county_1930.shp"
shp1    = r"" + path + "/tmp/CountiesND.shp"
shp2    = r"" + path + "/rawdata/NHGIS/CountiesND.shp"
shp3    = r"" + path + "/rawdata/NHGIS/CountiesND_WGS84.shp"
shp_nd  = r"" + path + "/rawdata/NHGIS/CountiesND.shp"
lyr1    = r"lyr1"
lyr2    = r"lyr2"

# Process: Make Feature Layer
print "I am creating a layer out of the original shapefile."
arcpy.MakeFeatureLayer_management(shp0, lyr1, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;DECADE DECADE VISIBLE NONE;NHGISNAM NHGISNAM VISIBLE NONE;NHGISST NHGISST VISIBLE NONE;NHGISCTY NHGISCTY VISIBLE NONE;ICPSRST ICPSRST VISIBLE NONE;ICPSRCTY ICPSRCTY VISIBLE NONE;ICPSRNAM ICPSRNAM VISIBLE NONE;STATENAM STATENAM VISIBLE NONE;ICPSRSTI ICPSRSTI VISIBLE NONE;ICPSRCTYI ICPSRCTYI VISIBLE NONE;ICPSRFIP ICPSRFIP VISIBLE NONE;STATE STATE VISIBLE NONE;COUNTY COUNTY VISIBLE NONE;PID PID VISIBLE NONE;X_CENTROID X_CENTROID VISIBLE NONE;Y_CENTROID Y_CENTROID VISIBLE NONE;GISJOIN GISJOIN VISIBLE NONE;GISJOIN2 GISJOIN2 VISIBLE NONE;SHAPE_AREA SHAPE_AREA VISIBLE NONE;SHAPE_LEN SHAPE_LEN VISIBLE NONE")

# Process: Add Join
print "I am joining the shapefile to the list of New Deal counties."
arcpy.AddJoin_management(lyr1, "GISJOIN", dbf, "GISJOIN", "KEEP_COMMON")

# Process: Dissolve
print "I am dissolving the shp & creating 1 obs for each New Deal county."
arcpy.Dissolve_management(lyr1, shp1, "Bridge.state;Bridge.stateicpsr;Bridge.countynd;Bridge.countyname", "", "MULTI_PART", "DISSOLVE_LINES")

# Process: Make Feature Layer
print "I am creating another layer out of the dissolved shapefile. This is needed to rename variable names."
arcpy.MakeFeatureLayer_management(shp1, lyr2, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;Bridge_sta state VISIBLE NONE;Bridge_s_1 stateicpsr VISIBLE NONE;Bridge_cou countynd VISIBLE NONE;Bridge_c_1 countynamend VISIBLE NONE")

# Process: Copy Features
print "I am exporting the layer as the final shapefile."
arcpy.CopyFeatures_management(lyr2, shp2, "", "0", "0", "0")

# Process: Add Geometry Attributes
print "I am adding area (in m2) and coordinates of the centroid."
arcpy.AddGeometryAttributes_management(shp2, "AREA", "", "SQUARE_METERS", "")
arcpy.AddGeometryAttributes_management(shp2, "CENTROID;CENTROID_INSIDE", "", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")

# Process: Project
print "I am projecting the map to WGS 1984."
arcpy.Project_management(shp2, shp3, "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]", "WGS_1984_(ITRF00)_To_NAD_1983", "PROJCS['USA_Contiguous_Albers_Equal_Area_Conic',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-96.0],PARAMETER['Standard_Parallel_1',29.5],PARAMETER['Standard_Parallel_2',45.5],PARAMETER['Latitude_Of_Origin',37.5],UNIT['Meter',1.0]]", "NO_PRESERVE_SHAPE", "")

# ---------------------------------------------------------------------------
#           2.      Compute intersections
#                 a. With New Deal counties
# ---------------------------------------------------------------------------
print "              Part 2: create the intersections"
print "                      a. New Deal counties"
for year in [1910,1940,1950]:
    print "Year " + str(year) + ":"
    print "   I am defining the variables."
    shp_10  = r"" + path + "/rawdata/NHGIS/US_county_" + str(year) + ".shp"
    shp_x   = r"" + path + "/tmp/intersection.shp"
    shp_ndx = r"" + path + "/rawdata/NHGIS/CountiesND_" + str(year) + ".shp"
    lyr     = r"lyr"

# Process: Intersect
    print "   I am intersecting the ND counties with the decadal counties."
    arcpy.Intersect_analysis(shp_nd + " #;" + shp_10 + " #", shp_x, "ALL", "", "INPUT")

# Process: Make Feature Layer
    print "   I am creating another layer out of the dissolved shapefile. This is needed to rename variable names."
    arcpy.MakeFeatureLayer_management(shp_x, lyr, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;FID_Counti FID_Counti HIDDEN NONE;state state VISIBLE NONE;stateicpsr stateicpsr VISIBLE NONE;countynd countynd VISIBLE NONE;countyname countyname VISIBLE NONE;POLY_AREA areand VISIBLE NONE;CENTROID_X CENTROID_X HIDDEN NONE;CENTROID_Y CENTROID_Y HIDDEN NONE;INSIDE_X INSIDE_X HIDDEN NONE;INSIDE_Y INSIDE_Y HIDDEN NONE;FID_US_cou FID_US_cou HIDDEN NONE;DECADE DECADE HIDDEN NONE;NHGISNAM NHGISNAM HIDDEN NONE;NHGISST NHGISST HIDDEN NONE;NHGISCTY NHGISCTY HIDDEN NONE;ICPSRST ICPSRST VISIBLE NONE;ICPSRCTY ICPSRCTY VISIBLE NONE;ICPSRNAM ICPSRNAM VISIBLE NONE;STATENAM STATENAM HIDDEN NONE;ICPSRSTI ICPSRSTI VISIBLE NONE;ICPSRCTYI ICPSRCTYI VISIBLE NONE;ICPSRFIP ICPSRFIP HIDDEN NONE;STATE_1 STATE_1 HIDDEN NONE;COUNTY COUNTY HIDDEN NONE;PID PID HIDDEN NONE;X_CENTROID X_CENTROID HIDDEN NONE;Y_CENTROID Y_CENTROID HIDDEN NONE;GISJOIN GISJOIN HIDDEN NONE;GISJOIN2 GISJOIN2 HIDDEN NONE;SHAPE_AREA SHAPE_AREA HIDDEN NONE;SHAPE_LEN SHAPE_LEN HIDDEN NONE")

# Process: Copy Features
    print "   I am exporting the layer as the final shapefile."
    arcpy.CopyFeatures_management(lyr, shp_ndx, "", "0", "0", "0")

# Process: Add Geometry Attributes
    print "   I am adding area (in m2)."
    arcpy.AddGeometryAttributes_management(shp_ndx, "AREA", "", "SQUARE_METERS", "")

# ---------------------------------------------------------------------------
#                 b. With Climatic divisions
# ---------------------------------------------------------------------------
print "                      b. Climatic divisions"
print "   I am defining the variables."
shp_nd  = r"" + path + "/rawdata/NHGIS/CountiesND.shp"
shp_40  = r"" + path + "/rawdata/NHGIS/US_county_1940.shp"
shp_cd  = r"" + path + "/rawdata/NOAA/OFFICIAL_CLIM_DIVISIONS.shp"
shp_x   = r"" + path + "/tmp/intersection.shp"
shp_ndx = r"" + path + "/rawdata/NOAA/CountiesND_CD.shp"
shp_40x = r"" + path + "/rawdata/NOAA/Counties1940_CD.shp"
lyr     = "lyr"

print "                          i. New Deal counties"
# Process: Intersect
print "   I am intersecting the ND counties with the climatic divisions."
arcpy.Intersect_analysis(shp_nd + " #;" + shp_cd + " #", shp_x, "ALL", "", "INPUT")

# Process: Make Feature Layer
print "   I am creating another layer out of the dissolved shapefile. This is needed to rename variables."
arcpy.MakeFeatureLayer_management(shp_x, lyr, "", "", "FID FID VISIBLE NONE;FID_Counti FID_Counti HIDDEN NONE;Shape Shape HIDDEN NONE;state state HIDDEN NONE;stateicpsr stateicpsr VISIBLE NONE;countynd countynd VISIBLE NONE;countyname countyname VISIBLE NONE;POLY_AREA POLY_AREA HIDDEN NONE;CENTROID_X CENTROID_X HIDDEN NONE;CENTROID_Y CENTROID_Y HIDDEN NONE;INSIDE_X INSIDE_X HIDDEN NONE;INSIDE_Y INSIDE_Y HIDDEN NONE;FID_OFFICI FID_OFFICI HIDDEN NONE;OBJECTID OBJECTID HIDDEN NONE;STATE_1 STATE_1 HIDDEN NONE;STATE_FIPS STATE_FIPS HIDDEN NONE;CD_2DIG CD_2DIG HIDDEN NONE;STATE_CODE STATE_CODE HIDDEN NONE;CLIMDIV CLIMDIV VISIBLE NONE;CD_NEW CD_NEW HIDDEN NONE;FIPS_CD FIPS_CD HIDDEN NONE;NCDC_GEO_I NCDC_GEO_I HIDDEN NONE;NAME NAME HIDDEN NONE;ST_ABBRV ST_ABBRV HIDDEN NONE;SHAPE_AREA SHAPE_AREA HIDDEN NONE;SHAPE_LEN SHAPE_LEN HIDDEN NONE;Area_mm Area_mm HIDDEN NONE")

# Process: Copy Features
print "   I am exporting the layer as the final shapefile."
arcpy.CopyFeatures_management(lyr, shp_ndx, "", "0", "0", "0")

# Process: Add Geometry Attributes
print "   I am adding area (in m2)."
arcpy.AddGeometryAttributes_management(shp_ndx, "AREA", "", "SQUARE_METERS", "")

print "                         ii. 1940 counties"
# Process: Intersect
print "   I am intersecting the 1940 counties with the climatic divisions."
arcpy.Intersect_analysis(shp_40 + " #;" + shp_cd + " #", shp_x, "ALL", "", "INPUT")

# Process: Make Feature Layer
print "   I am creating another layer out of the dissolved shapefile. This is needed to rename variables."
arcpy.MakeFeatureLayer_management(shp_x, lyr, "", "", "FID FID VISIBLE NONE;FID_Counti FID_Counti HIDDEN NONE;Shape Shape HIDDEN NONE;state state HIDDEN NONE;stateicpsr stateicpsr VISIBLE NONE;countynd countynd VISIBLE NONE;countyname countyname VISIBLE NONE;POLY_AREA POLY_AREA HIDDEN NONE;CENTROID_X CENTROID_X HIDDEN NONE;CENTROID_Y CENTROID_Y HIDDEN NONE;INSIDE_X INSIDE_X HIDDEN NONE;INSIDE_Y INSIDE_Y HIDDEN NONE;FID_OFFICI FID_OFFICI HIDDEN NONE;OBJECTID OBJECTID HIDDEN NONE;STATE_1 STATE_1 HIDDEN NONE;STATE_FIPS STATE_FIPS HIDDEN NONE;CD_2DIG CD_2DIG HIDDEN NONE;STATE_CODE STATE_CODE HIDDEN NONE;CLIMDIV CLIMDIV VISIBLE NONE;CD_NEW CD_NEW HIDDEN NONE;FIPS_CD FIPS_CD HIDDEN NONE;NCDC_GEO_I NCDC_GEO_I HIDDEN NONE;NAME NAME HIDDEN NONE;ST_ABBRV ST_ABBRV HIDDEN NONE;SHAPE_AREA SHAPE_AREA HIDDEN NONE;SHAPE_LEN SHAPE_LEN HIDDEN NONE;Area_mm Area_mm HIDDEN NONE")

# Process: Copy Features
print "   I am exporting the layer as the final shapefile."
arcpy.CopyFeatures_management(lyr, shp_40x, "", "0", "0", "0")

# Process: Add Geometry Attributes
print "   I am adding area (in m2)."
arcpy.AddGeometryAttributes_management(shp_40x, "AREA", "", "SQUARE_METERS", "")

# ---------------------------------------------------------------------------
#                 c. With Electoral districts (Committees)
# ---------------------------------------------------------------------------

# Local variables:

dist62   = r"" + path + "/rawdata/Lewis/districts062.shp"
dist73   = r"" + path + "/rawdata/Lewis/districts073.shp"
shp_nd   = r"" + path + "/rawdata/NHGIS/CountiesND.shp"
shp_x    = r"" + path + "/tmp/intersection.shp"
shp_nd62 = r"" + path + "/rawdata/Lewis/CountiesND_62.shp"
shp_nd73 = r"" + path + "/rawdata/Lewis/CountiesND_73.shp"

lyr     = r"lyr"


print "I am creating a correspondence between counties and 62th Congress"
# Process: Intersect
arcpy.Intersect_analysis(dist62 + " #;" + shp_nd + " #", shp_x, "ALL", "", "INPUT")

# Process: Make Feature Layer
print "   I am creating another layer out of the dissolved shapefile. This is needed to rename variable names."
arcpy.MakeFeatureLayer_management(shp_x, lyr, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;FID_Counti FID_Counti HIDDEN NONE;state state VISIBLE NONE;stateicpsr stateicpsr VISIBLE NONE;countynd countynd VISIBLE NONE;countyname countyname VISIBLE NONE;POLY_AREA areand VISIBLE NONE;CENTROID_X CENTROID_X HIDDEN NONE;CENTROID_Y CENTROID_Y HIDDEN NONE;INSIDE_X INSIDE_X HIDDEN NONE;INSIDE_Y INSIDE_Y HIDDEN NONE;FID_US_cou FID_US_cou HIDDEN NONE;DECADE DECADE HIDDEN NONE;NHGISNAM NHGISNAM HIDDEN NONE;NHGISST NHGISST HIDDEN NONE;NHGISCTY NHGISCTY HIDDEN NONE;ICPSRST ICPSRST VISIBLE NONE;ICPSRCTY ICPSRCTY VISIBLE NONE;ICPSRNAM ICPSRNAM VISIBLE NONE;STATENAM STATENAM HIDDEN NONE;ICPSRSTI ICPSRSTI VISIBLE NONE;ICPSRCTYI ICPSRCTYI VISIBLE NONE;ICPSRFIP ICPSRFIP HIDDEN NONE;STATE_1 STATE_1 HIDDEN NONE;COUNTY COUNTY HIDDEN NONE;PID PID HIDDEN NONE;X_CENTROID X_CENTROID HIDDEN NONE;Y_CENTROID Y_CENTROID HIDDEN NONE;GISJOIN GISJOIN HIDDEN NONE;GISJOIN2 GISJOIN2 HIDDEN NONE;SHAPE_AREA SHAPE_AREA HIDDEN NONE;SHAPE_LEN SHAPE_LEN HIDDEN NONE")

# Process: Copy Features
print "   I am exporting the layer as the final shapefile."
arcpy.CopyFeatures_management(lyr, shp_nd62, "", "0", "0", "0")

# Process: Add Geometry Attributes
print "   I am adding area (in m2)."
arcpy.AddGeometryAttributes_management(shp_nd62, "AREA_GEODESIC", "", "SQUARE_METERS", "")


print "I am creating a correspondence between counties and 73rd Congress"
# Process: Intersect
arcpy.Intersect_analysis(dist73 + " #;" + shp_nd + " #", shp_x, "ALL", "", "INPUT")

# Process: Make Feature Layer
print "   I am creating another layer out of the dissolved shapefile. This is needed to rename variable names."
arcpy.MakeFeatureLayer_management(shp_x, lyr, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;FID_Counti FID_Counti HIDDEN NONE;state state VISIBLE NONE;stateicpsr stateicpsr VISIBLE NONE;countynd countynd VISIBLE NONE;countyname countyname VISIBLE NONE;POLY_AREA areand VISIBLE NONE;CENTROID_X CENTROID_X HIDDEN NONE;CENTROID_Y CENTROID_Y HIDDEN NONE;INSIDE_X INSIDE_X HIDDEN NONE;INSIDE_Y INSIDE_Y HIDDEN NONE;FID_US_cou FID_US_cou HIDDEN NONE;DECADE DECADE HIDDEN NONE;NHGISNAM NHGISNAM HIDDEN NONE;NHGISST NHGISST HIDDEN NONE;NHGISCTY NHGISCTY HIDDEN NONE;ICPSRST ICPSRST VISIBLE NONE;ICPSRCTY ICPSRCTY VISIBLE NONE;ICPSRNAM ICPSRNAM VISIBLE NONE;STATENAM STATENAM HIDDEN NONE;ICPSRSTI ICPSRSTI VISIBLE NONE;ICPSRCTYI ICPSRCTYI VISIBLE NONE;ICPSRFIP ICPSRFIP HIDDEN NONE;STATE_1 STATE_1 HIDDEN NONE;COUNTY COUNTY HIDDEN NONE;PID PID HIDDEN NONE;X_CENTROID X_CENTROID HIDDEN NONE;Y_CENTROID Y_CENTROID HIDDEN NONE;GISJOIN GISJOIN HIDDEN NONE;GISJOIN2 GISJOIN2 HIDDEN NONE;SHAPE_AREA SHAPE_AREA HIDDEN NONE;SHAPE_LEN SHAPE_LEN HIDDEN NONE")

# Process: Copy Features
print "   I am exporting the layer as the final shapefile."
arcpy.CopyFeatures_management(lyr, shp_nd73, "", "0", "0", "0")

# Process: Add Geometry Attributes
print "   I am adding area (in m2)."
arcpy.AddGeometryAttributes_management(shp_nd73, "AREA_GEODESIC", "", "SQUARE_METERS", "")


# ---------------------------------------------------------------------------
#           3.      Erase junk
# ---------------------------------------------------------------------------
extensions = ["cpg","dbf","prj","sbn","sbx","shp","shp.xml","shx"]
files      = ["intersection","CountiesND"]

print "I am erasing the junk I have created."
for file in files:
    for extension in extensions:
        os.remove(r"" + path + "/tmp/" + file + "." + extension)

try:
    os.rmdir(r"" + path + "/tmp/")
except:
    pass
